diff --git a/.gitignore b/.gitignore index f43115e..e1ba14f 100644 --- a/.gitignore +++ b/.gitignore @@ -83,3 +83,4 @@ venv/ **/_version.py .conda/ +debug/* diff --git a/derotation/analysis/fit_ellipse.py b/derotation/analysis/fit_ellipse.py new file mode 100644 index 0000000..f92a887 --- /dev/null +++ b/derotation/analysis/fit_ellipse.py @@ -0,0 +1,165 @@ +import logging +from pathlib import Path +from typing import Tuple + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.patches import Ellipse +from scipy.optimize import OptimizeResult, least_squares + + +def fit_ellipse_to_points( + centers: np.ndarray, +) -> Tuple[int, int, int, int, int]: + """Fit an ellipse to the points using least squares optimization. + + Parameters + ---------- + centers : np.ndarray + The centers of the largest blob in each image. + + Returns + ------- + Tuple[int, int, int, int, int] + The center of the ellipse (center_x, center_y), the semi-major + axis (a), the semi-minor axis (b), and the rotation angle (theta). + """ + # Convert centers to numpy array + centers = np.array(centers) + x = centers[:, 0] + y = centers[:, 1] + + # Find the extreme points for the initial ellipse estimate + topmost = centers[np.argmin(y)] + rightmost = centers[np.argmax(x)] + bottommost = centers[np.argmax(y)] + leftmost = centers[np.argmin(x)] + + # Initial parameters: (center_x, center_y, semi_major_axis, + # semi_minor_axis, rotation_angle) + initial_center = np.mean( + [topmost, bottommost, leftmost, rightmost], axis=0 + ) + semi_major_axis = np.linalg.norm(rightmost - leftmost) / 2 + semi_minor_axis = np.linalg.norm(topmost - bottommost) / 2 + rotation_angle = 0 # Start with no rotation + + initial_params = [ + initial_center[0], + initial_center[1], + semi_major_axis, + semi_minor_axis, + rotation_angle, + ] + + logging.info("Fitting ellipse to points...") + + # Objective function to minimize: sum of squared distances to ellipse + def ellipse_residuals(params, x, y): + center_x, center_y, a, b, theta = params # theta is in radians + cos_angle = np.cos(theta) + sin_angle = np.sin(theta) + + # Rotate the points to align with the ellipse axes + x_rot = cos_angle * (x - center_x) + sin_angle * (y - center_y) + y_rot = -sin_angle * (x - center_x) + cos_angle * (y - center_y) + + # Ellipse equation: (x_rot^2 / a^2) + (y_rot^2 / b^2) = 1 + return (x_rot / a) ** 2 + (y_rot / b) ** 2 - 1 + + # Use least squares optimization to fit the ellipse to the points + result: OptimizeResult = least_squares( + ellipse_residuals, + initial_params, + args=(x, y), + loss="huber", # minimize the influence of outliers + ) + + # Extract optimized parameters + center_x, center_y, a, b, theta = result.x + + return center_x, center_y, a, b, theta + + +def plot_ellipse_fit_and_centers( + centers: np.ndarray, + center_x: int, + center_y: int, + a: int, + b: int, + theta: int, + image_stack: np.ndarray, + debug_plots_folder: Path, + saving_name: str = "ellipse_fit.png", +): + """Plot the fitted ellipse on the largest blob centers. + + Parameters + ---------- + centers : np.ndarray + The centers of the largest blob in each image. + center_x : int + The x-coordinate of the center of the ellipse + center_y : int + The y-coordinate of the center of the ellipse + a : int + The semi-major axis of the ellipse + b : int + 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) + x = centers[:, 0] + y = centers[:, 1] + + # Create the plot + fig, ax = plt.subplots(figsize=(8, 8)) + + max_projection = image_stack.max(axis=0) + # plot behind a frame of the original image + ax.imshow(max_projection, cmap="gray") + ax.scatter(x, y, label="Largest Blob Centers", color="red") + + # Plot fitted ellipse + ellipse = Ellipse( + (center_x, center_y), + width=2 * a, + height=2 * b, + angle=np.degrees(theta), + edgecolor="blue", + facecolor="none", + label="Fitted Ellipse", + ) + ax.add_patch(ellipse) + + # Plot center of fitted ellipse + ax.scatter( + center_x, + center_y, + color="green", + marker="x", + s=100, + label="Ellipse Center", + ) + + # Add some plot formatting + ax.set_xlim(0, image_stack.shape[1]) + ax.set_ylim( + image_stack.shape[1], 0 + ) # Invert y-axis to match image coordinate system + ax.set_aspect("equal") + ax.legend() + ax.grid(True) + ax.set_title("Fitted Ellipse on largest blob centers") + ax.axis("off") + + # save + plt.savefig(debug_plots_folder / saving_name) diff --git a/derotation/analysis/incremental_derotation_pipeline.py b/derotation/analysis/incremental_derotation_pipeline.py index 25c1e3e..f0f8873 100644 --- a/derotation/analysis/incremental_derotation_pipeline.py +++ b/derotation/analysis/incremental_derotation_pipeline.py @@ -6,12 +6,14 @@ import matplotlib.pyplot as plt import numpy as np import pandas as pd -from matplotlib.patches import Ellipse from scipy.ndimage import rotate -from scipy.optimize import OptimizeResult, least_squares from skimage.feature import blob_log from tqdm import tqdm +from derotation.analysis.fit_ellipse import ( + fit_ellipse_to_points, + plot_ellipse_fit_and_centers, +) from derotation.analysis.full_derotation_pipeline import FullPipeline @@ -489,10 +491,23 @@ def find_center_of_rotation(self) -> Tuple[int, int]: ) # Fit an ellipse to the largest blob centers and get its center - center_x, center_y, a, b, theta = self.fit_ellipse_to_points( + center_x, center_y, a, b, theta = fit_ellipse_to_points( coord_first_blob_of_every_image ) + if self.debugging_plots: + plot_ellipse_fit_and_centers( + coord_first_blob_of_every_image, + center_x, + center_y, + a, + b, + theta, + image_stack=mean_images, + debug_plots_folder=self.debug_plots_folder, + saving_name="ellipse_fit.png", + ) + logging.info( f"Center of ellipse: ({center_x:.2f}, {center_y:.2f}), " f"semi-major axis: {a:.2f}, semi-minor axis: {b:.2f}" @@ -531,9 +546,15 @@ def get_coords_of_largest_blob( blobs[i][blobs[i][:, 2].argsort()] for i in range(len(image_stack)) ] - coord_first_blob_of_every_image = [ - blobs[i][0][:2].astype(int) for i in range(len(image_stack)) - ] + try: + coord_first_blob_of_every_image = [ + blobs[i][0][:2].astype(int) for i in range(len(image_stack)) + ] + except IndexError: + raise ValueError( + "No blobs were found. Try changing the parameters of the " + + "blob detection algorithm." + ) # invert x, y order coord_first_blob_of_every_image = [ @@ -573,150 +594,3 @@ def plot_blob_detection(self, blobs: list, image_stack: np.ndarray): # save the plot plt.savefig(self.debug_plots_folder / "blobs.png") - - def fit_ellipse_to_points( - self, centers: np.ndarray - ) -> Tuple[int, int, int, int, int]: - """Fit an ellipse to the points using least squares optimization. - - Parameters - ---------- - centers : np.ndarray - The centers of the largest blob in each image. - - Returns - ------- - Tuple[int, int, int, int, int] - The center of the ellipse (center_x, center_y), the semi-major - axis (a), the semi-minor axis (b), and the rotation angle (theta). - """ - # Convert centers to numpy array - centers = np.array(centers) - x = centers[:, 0] - y = centers[:, 1] - - # Find the extreme points for the initial ellipse estimate - topmost = centers[np.argmin(y)] - rightmost = centers[np.argmax(x)] - bottommost = centers[np.argmax(y)] - leftmost = centers[np.argmin(x)] - - # Initial parameters: (center_x, center_y, semi_major_axis, - # semi_minor_axis, rotation_angle) - initial_center = np.mean( - [topmost, bottommost, leftmost, rightmost], axis=0 - ) - semi_major_axis = np.linalg.norm(rightmost - leftmost) / 2 - semi_minor_axis = np.linalg.norm(topmost - bottommost) / 2 - rotation_angle = 0 # Start with no rotation - - initial_params = [ - initial_center[0], - initial_center[1], - semi_major_axis, - semi_minor_axis, - rotation_angle, - ] - - logging.info("Fitting ellipse to points...") - - # Objective function to minimize: sum of squared distances to ellipse - def ellipse_residuals(params, x, y): - center_x, center_y, a, b, theta = params - cos_angle = np.cos(theta) - sin_angle = np.sin(theta) - - # Rotate the points to align with the ellipse axes - x_rot = cos_angle * (x - center_x) + sin_angle * (y - center_y) - y_rot = -sin_angle * (x - center_x) + cos_angle * (y - center_y) - - # Ellipse equation: (x_rot^2 / a^2) + (y_rot^2 / b^2) = 1 - return (x_rot / a) ** 2 + (y_rot / b) ** 2 - 1 - - # Use least squares optimization to fit the ellipse to the points - result: OptimizeResult = least_squares( - ellipse_residuals, initial_params, args=(x, y) - ) - - # Extract optimized parameters - center_x, center_y, a, b, theta = result.x - - if self.debugging_plots: - self.plot_ellipse_fit_and_centers( - centers, center_x, center_y, a, b, theta - ) - - return center_x, center_y, a, b, theta - - def plot_ellipse_fit_and_centers( - self, - centers: np.ndarray, - center_x: int, - center_y: int, - a: int, - b: int, - theta: int, - ): - """Plot the fitted ellipse on the largest blob centers. - - Parameters - ---------- - centers : np.ndarray - The centers of the largest blob in each image. - center_x : int - The x-coordinate of the center of the ellipse - center_y : int - The y-coordinate of the center of the ellipse - a : int - The semi-major axis of the ellipse - b : int - The semi-minor axis of the ellipse - theta : int - The rotation angle of the ellipse - """ - # Convert centers to numpy array - centers = np.array(centers) - x = centers[:, 0] - y = centers[:, 1] - - # Create the plot - fig, ax = plt.subplots(figsize=(8, 8)) - # plot behind a frame of the original image - ax.imshow(self.image_stack[0], cmap="gray") - ax.scatter(x, y, label="Largest Blob Centers", color="red") - - # Plot fitted ellipse - ellipse = Ellipse( - (center_x, center_y), - width=2 * a, - height=2 * b, - angle=np.degrees(theta), - edgecolor="blue", - facecolor="none", - label="Fitted Ellipse", - ) - ax.add_patch(ellipse) - - # Plot center of fitted ellipse - ax.scatter( - center_x, - center_y, - color="green", - marker="x", - s=100, - label="Ellipse Center", - ) - - # Add some plot formatting - ax.set_xlim(0, self.image_stack.shape[1]) - ax.set_ylim( - self.image_stack.shape[1], 0 - ) # Invert y-axis to match image coordinate system - ax.set_aspect("equal") - ax.legend() - ax.grid(True) - ax.set_title("Fitted Ellipse on largest blob centers") - ax.axis("off") - - # save - plt.savefig(self.debug_plots_folder / "ellipse_fit.png") diff --git a/derotation/simulate/basic_rotator.py b/derotation/simulate/basic_rotator.py deleted file mode 100644 index 26ace59..0000000 --- a/derotation/simulate/basic_rotator.py +++ /dev/null @@ -1,153 +0,0 @@ -import copy -from typing import Optional, Tuple - -import numpy as np -from scipy.ndimage import affine_transform - - -class Rotator: - def __init__( - self, - angles: np.ndarray, - image_stack: np.ndarray, - center: Optional[Tuple[int, int]] = None, - ) -> 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] - - if center is None: - self.center = np.array(image_stack.shape[1:]) / 2 - else: - self.center = np.array(center) - - 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) - - 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) - rotated_image_stack[i][j] = rotated_frame[j] - - return rotated_image_stack - - def rotate(self, image: np.ndarray, angle: 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) - - # Rotation matrix clockwise if angle is positive - rotation_matrix = np.array( - [ - [cos, -sin], - [sin, cos], - ] - ) - - # Compute offset so rotation is around the center - offset = self.center - rotation_matrix @ self.center - - # 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=self.get_blank_pixels_value(), - ) - - 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) diff --git a/derotation/simulate/line_scanning_microscope.py b/derotation/simulate/line_scanning_microscope.py new file mode 100644 index 0000000..3298391 --- /dev/null +++ b/derotation/simulate/line_scanning_microscope.py @@ -0,0 +1,318 @@ +from typing import Optional, Tuple + +import numpy as np +import tqdm +from scipy.ndimage import affine_transform + + +class Rotator: + def __init__( + self, + angles: np.ndarray, + image_stack: np.ndarray, + center: Optional[Tuple[int, int]] = 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 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. + + 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 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. + + 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. + 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 : float, optional + The z angle of the rotation plane in degrees in relation to the + scanning plane. If 0, the rotation plane is the same as + the scanning plane, by default None. + rotation_plane_orientation : float, optional + The angle of the rotation plane in the x-y plane in degrees, + 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) + 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], ( + 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]})" + ) + + 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." + ) + + self.image_stack = image_stack + self.num_lines_per_frame = image_stack.shape[1] + + if center is None: + self.center = np.array(image_stack.shape[1:]) / 2 + else: + self.center = np.array(center) + + 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.create_homography_matrices() + self.calculate_pixel_shift() + print(f"Pixel shift: {self.ps}") + print(f"New image size: {self.image_size}") + + # reshape the angles to the shape of the image + # stack to ease indexing + self.angles = angles[: image_stack.shape[0] * self.image_size].reshape( + image_stack.shape[0], self.image_size + ) + print("New angles shape:", self.angles.shape) + + if blank_pixel_val is None: + self.blank_pixel_val = self.get_blank_pixels_value() + else: + 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( + [ + [1, 0, 0], + [0, np.cos(self.rotation_plane_angle), 0], + [0, 0, 1], + ] + ) + + # 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 = line_length - np.round( + np.abs(line_length * np.cos(self.rotation_plane_angle)) + ).astype(int) + + # round to the nearest even number + self.ps += self.ps % 2 + + # final image size depends on the pixel shift + 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: + return image[ + # centered in the rows + self.ps // 2 : -self.ps // 2, + # take the left side of the image + : self.image_size, + ] + + 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 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. + + 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 = np.empty( + (self.image_stack.shape[0], self.image_size, self.image_size), + dtype=np.float64, + ) + + for i, image in tqdm.tqdm( + enumerate(self.image_stack), total=self.image_stack.shape[0] + ): + if np.all( + # don't bother if rotation is less than 0.01 degrees + np.isclose(self.angles[i], 0, atol=1e-2) + ): + 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) + + # store the rotated image line + rotated_image_stack[i][j] = rotated_image[j] + + return rotated_image_stack + + def homography_rotation_to_scanning_plane( + self, image: np.ndarray + ) -> np.ndarray: + """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. + + Returns + ------- + np.ndarray + The transformed image. + """ + + # 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) + + return image + + def rotate_sample(self, image: np.ndarray, angle: 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) + + # Rotation matrix clockwise if angle is positive + rotation_matrix = np.array( + [ + [cos, -sin], + [sin, cos], + ] + ) + + # Compute offset so rotation is around the center + offset = self.center - rotation_matrix @ self.center + + # 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=self.blank_pixel_val, + ) + + 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) diff --git a/examples/elliptical_rotations.py b/examples/elliptical_rotations.py new file mode 100644 index 0000000..69ab8ae --- /dev/null +++ b/examples/elliptical_rotations.py @@ -0,0 +1,104 @@ +# 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 + +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") + + +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 +) +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/examples/rotate_and_derotate_a_square.py b/examples/rotate_and_derotate_a_square.py index dc7ec52..36834b5 100644 --- a/examples/rotate_and_derotate_a_square.py +++ b/examples/rotate_and_derotate_a_square.py @@ -2,7 +2,7 @@ import numpy as np from derotation.derotate_by_line import derotate_an_image_array_line_by_line -from derotation.simulate.basic_rotator import Rotator +from derotation.simulate.line_scanning_microscope import Rotator # make a simple image, a square in a black background image = np.empty((100, 100)) diff --git a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py b/tests/test_integration/test_finding_center_of_rotation.py similarity index 95% rename from tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py rename to tests/test_integration/test_finding_center_of_rotation.py index a5703f9..4d229ef 100644 --- a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py +++ b/tests/test_integration/test_finding_center_of_rotation.py @@ -53,7 +53,7 @@ IncrementalPipeline, ) from derotation.derotate_by_line import derotate_an_image_array_line_by_line -from derotation.simulate.basic_rotator import Rotator +from derotation.simulate.line_scanning_microscope import Rotator # ----------------------------------------------------- # Prepare the 3D image stack and the rotation angles @@ -63,8 +63,11 @@ def create_sample_image_with_two_cells( center_of_bright_cell: Tuple[int, int] = (50, 10), center_of_dimmer_cell: Tuple[int, int] = (60, 60), + lines_per_frame: int = 100, + 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 @@ -76,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 ------- @@ -84,16 +93,11 @@ def create_sample_image_with_two_cells( """ # Initialize a black image of size 100x100 - image = np.zeros((100, 100), dtype=np.uint8) + image = np.zeros((lines_per_frame, lines_per_frame), dtype=np.uint8) # Define the circle's parameters - radius = 5 # radius of the circle white_value = 255 # white color for the circle - # add an extra gray circle at the bottom right - radius2 = 5 - gray_value = 128 - # Draw a white circle in the top center y, x = np.ogrid[: image.shape[0], : image.shape[1]] mask = (x - center_of_bright_cell[0]) ** 2 + ( @@ -101,11 +105,14 @@ def create_sample_image_with_two_cells( ) ** 2 <= radius**2 image[mask] = white_value - # Draw a gray circle in the bottom right - mask2 = (x - center_of_dimmer_cell[0]) ** 2 + ( - y - center_of_dimmer_cell[1] - ) ** 2 <= radius2**2 - image[mask2] = gray_value + if second_cell: + # add an extra gray circle at the bottom right + gray_value = 128 + # Draw a gray circle in the bottom right + mask2 = (x - center_of_dimmer_cell[0]) ** 2 + ( + y - center_of_dimmer_cell[1] + ) ** 2 <= radius**2 + image[mask2] = gray_value return image @@ -471,7 +478,9 @@ def test_blob_detection_on_derotated_stack( # If not, the derotation was not successful and the test fails # Detect the blobs in the derotated stack - blobs = [blob_log(img) for img in derotated_sinusoidal] + blobs = [ + blob_log(img, min_sigma=3, max_sigma=5) for img in derotated_sinusoidal + ] # Get the center of the blobs # for every frame, place first the blob with the smallest x value diff --git a/tests/test_integration/test_rotation_out_of_plane.py b/tests/test_integration/test_rotation_out_of_plane.py new file mode 100644 index 0000000..2a84189 --- /dev/null +++ b/tests/test_integration/test_rotation_out_of_plane.py @@ -0,0 +1,178 @@ +from pathlib import Path +from typing import Tuple + +import numpy as np +import pytest + +from derotation.analysis.fit_ellipse import ( + fit_ellipse_to_points, + plot_ellipse_fit_and_centers, +) +from derotation.simulate.line_scanning_microscope import Rotator +from tests.test_integration.test_finding_center_of_rotation import ( + create_image_stack, + create_rotation_angles, + create_sample_image_with_two_cells, +) + + +def rotate_image_stack( + plane_angle: float = 0, + num_frames: int = 50, + pad: int = 20, + 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 : 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 : float, optional + Orientation of the rotation plane in degrees, by default None. + + Returns + ------- + Tuple[np.ndarray, np.ndarray, Rotator, int] + Original image stack, rotated image stack, Rotator instance, + and number of frames. + """ + cells = create_sample_image_with_two_cells( + lines_per_frame=100, second_cell=False, radius=1 + ) + cells = np.pad(cells, ((pad, pad), (pad, pad)), mode="constant") + cells[cells == 0] = 80 + + image_stack = create_image_stack(cells, num_frames=num_frames) + _, angles = create_rotation_angles(image_stack.shape) + + rotator = Rotator( + angles, + image_stack, + rotation_plane_angle=plane_angle, + blank_pixel_val=0, + rotation_plane_orientation=orientation, + ) + rotated_image_stack = rotator.rotate_by_line() + return image_stack, rotated_image_stack, rotator, num_frames + + +@pytest.mark.parametrize( + "plane_angle,exp_orientation", + [ + (0, 0), + (15, 0), + (25, 0), + (25, 45), + (25, 90), + (40, 0), + (40, 45), + (40, 90), + ], +) +def test_max_projection( + plane_angle: float, + exp_orientation: float, + atol: int = 15, +) -> None: + """ + Tests if the max projection of an image stack fits an ellipse, + and verify orientation and major/minor axes are close to expected values: + - if plane_angle is 0, the major and minor axes should be close to each + other; + - if exp_orientation is not None, the orientation should be close to + it. + - if exp_orientation is None, the orientation should be close to 0. + + Giving by default an allowed tolerance of 15 degrees as I care about + the test tracking the general shape of the ellipse, not the exact values. + + Parameters + ---------- + plane_angle : float + The angle of the rotation plane in degrees. + exp_orientation : float + Expected orientation of the rotation plane, by default None. + atol : int, optional + Allowed tolerance for orientation, by default 15. + """ + _, rotated_image_stack, *_ = rotate_image_stack( + plane_angle=plane_angle, orientation=exp_orientation + ) + + # get the location of the brightest pixel in any frame as a center + centers = np.array( + [ + np.unravel_index(np.argmax(frame), frame.shape) + for frame in rotated_image_stack + ] + ) + # invert centers to (x, y) format + centers = np.array([(x, y) for y, x in centers]) + + xc, yc, a, b, orientation = fit_ellipse_to_points(centers) + + if exp_orientation is None and plane_angle == 0: + # lower tolerance for the major and minor axes + if not np.allclose(a, b, atol=atol - 5): + plot_ellipse_fit_and_centers( + image_stack=rotated_image_stack, + centers=centers, + center_x=xc, + center_y=yc, + a=a, + b=b, + theta=orientation, + debug_plots_folder=Path("debug/"), + saving_name=f"ellipse_fit_{plane_angle}_{exp_orientation}.png", + ) + assert ( + False + ), f"Major and minor axes should be close, instead got {a} and {b}" + elif exp_orientation is not None: + # Major axis orientation in clockwise direction as radians. + exp_orientation = np.deg2rad(exp_orientation) + if not np.allclose(orientation, exp_orientation, atol=atol): + plot_ellipse_fit_and_centers( + image_stack=rotated_image_stack, + centers=centers, + center_x=xc, + center_y=yc, + a=a, + b=b, + theta=orientation, + debug_plots_folder=Path("debug/"), + saving_name=f"ellipse_fit_{plane_angle}_{exp_orientation}.png", + ) + assert False, ( + f"Orientation should be close to {exp_orientation}, " + f"instead got {orientation}" + ) + else: + # check that the major axis (a) differs from the minor axis (b) + # as calculated with cosine + decrease = np.cos(np.deg2rad(plane_angle)) + expected_b = a - a * decrease + # lower tolerance for the major and minor axes + if np.allclose(b, expected_b, atol=atol - 5): + plot_ellipse_fit_and_centers( + image_stack=rotated_image_stack, + centers=centers, + center_x=xc, + center_y=yc, + a=a, + b=b, + theta=orientation, + debug_plots_folder=Path("debug/"), + saving_name=f"ellipse_fit_{plane_angle}_{exp_orientation}.png", + ) + assert False, ( + f"Difference between major and minor axes should be close to " + f"{expected_b}, instead got {b}" + ) diff --git a/tests/test_regression/recreate_target/shared.py b/tests/test_regression/recreate_target/shared.py index e7ef854..d4430a8 100644 --- a/tests/test_regression/recreate_target/shared.py +++ b/tests/test_regression/recreate_target/shared.py @@ -5,7 +5,7 @@ from PIL import Image from derotation.derotate_by_line import derotate_an_image_array_line_by_line -from derotation.simulate.basic_rotator import Rotator +from derotation.simulate.line_scanning_microscope import Rotator NUMBER_OF_FRAMES = 3 ROTATED_IMAGES_PATH = Path("tests/test_regression/images/rotator") diff --git a/tests/test_regression/test_basic_rotator_with_different_center_of_rotation.py b/tests/test_regression/test_basic_rotator_with_different_center_of_rotation.py index 2b7f5f6..e8ca73a 100644 --- a/tests/test_regression/test_basic_rotator_with_different_center_of_rotation.py +++ b/tests/test_regression/test_basic_rotator_with_different_center_of_rotation.py @@ -3,7 +3,7 @@ from assertions import compare_images from PIL import Image -from derotation.simulate.basic_rotator import Rotator +from derotation.simulate.line_scanning_microscope import Rotator from tests.test_regression.recreate_target.shared import ( ROTATED_IMAGES_PATH, center_formatting, diff --git a/tests/test_unit/test_rotator.py b/tests/test_unit/test_rotator.py index 4879e97..6f5f822 100644 --- a/tests/test_unit/test_rotator.py +++ b/tests/test_unit/test_rotator.py @@ -1,6 +1,6 @@ import numpy as np -from derotation.simulate.basic_rotator import Rotator +from derotation.simulate.line_scanning_microscope import Rotator def test_Rotator_constructor():