From f301d96a1da91a3f7d859ae0a4bea8dfc97a7b12 Mon Sep 17 00:00:00 2001 From: lauraporta <29216006+lauraporta@users.noreply.github.com> Date: Thu, 5 Dec 2024 17:57:16 +0000 Subject: [PATCH] =?UTF-8?q?WIP=20tests=20-=20breaking=20=F0=9F=90=9B?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ...plane.py => test_rotation_out_of_plane.py} | 147 +++++++++++++++--- 1 file changed, 127 insertions(+), 20 deletions(-) rename tests/test_integration/{rotation_out_of_plane.py => test_rotation_out_of_plane.py} (51%) diff --git a/tests/test_integration/rotation_out_of_plane.py b/tests/test_integration/test_rotation_out_of_plane.py similarity index 51% rename from tests/test_integration/rotation_out_of_plane.py rename to tests/test_integration/test_rotation_out_of_plane.py index d9e68b0..e951cf5 100644 --- a/tests/test_integration/rotation_out_of_plane.py +++ b/tests/test_integration/test_rotation_out_of_plane.py @@ -1,22 +1,24 @@ from typing import Optional, Tuple import numpy as np +import pytest from matplotlib import pyplot as plt +from skimage.draw import ellipse_perimeter from skimage.feature import canny from skimage.transform import hough_ellipse -from test_finding_center_of_rotation_by_joining_two_pipelines import ( + +from derotation.simulate.line_scanning_microscope import Rotator +from tests.test_integration.test_finding_center_of_rotation_by_joining_two_pipelines import ( create_image_stack, create_rotation_angles, create_sample_image_with_two_cells, ) -from derotation.simulate.line_scanning_microscope import Rotator - def rotate_image_stack( - plane_angle: int, - num_frames: int, - pad: int, + plane_angle: int = 0, + num_frames: int = 50, + pad: int = 20, orientation: Optional[int] = None, ) -> Tuple[np.ndarray, np.ndarray, Rotator, int]: """ @@ -125,42 +127,147 @@ def make_plot( plt.savefig(f"debug/{title}.png") +@pytest.mark.parametrize( + "plane_angle,expected_orientation", + [ + (0, None), + (15, None), + (25, None), + (25, 45), + (25, 90), + (25, 135), + (40, None), + (40, 45), + (40, 90), + (40, 135), + ], +) def test_max_projection( - rotated_image_stack: np.ndarray, - expected_orientation: Optional[int] = None, - atol: int = 5, + plane_angle: int, + expected_orientation: Optional[int], + atol: int = 15, ) -> None: """ - Tests the max projection for an image stack to validate rotation. + 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 expected_orientation is not None, the orientation should be close to it. + - if expected_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 ---------- - rotated_image_stack : np.ndarray - Rotated image stack. + plane_angle : int + The angle of the rotation plane in degrees. expected_orientation : int, optional Expected orientation of the rotation plane, by default None. atol : int, optional Allowed tolerance for orientation, by default 5. """ + _, rotated_image_stack, *_ = rotate_image_stack( + plane_angle=plane_angle, orientation=expected_orientation + ) + max_projection = rotated_image_stack.max(axis=0) - edges = canny(max_projection, sigma=7, low_threshold=35, high_threshold=50) + edges = canny(max_projection, sigma=7, low_threshold=30, high_threshold=50) result = hough_ellipse( - edges, accuracy=10, threshold=100, min_size=30, max_size=100 + edges, accuracy=9, threshold=110, min_size=35, max_size=100 ) result.sort(order="accumulator") best = result[-1] - yc, xc, a, b, o = [int(best[i]) for i in range(1, len(best))] + yc, xc, a, b, orientation = [int(best[i]) for i in range(1, len(best))] - if expected_orientation is not None: - assert np.allclose( - o, expected_orientation, atol=atol - ), f"Orientation should be close to {expected_orientation}" + # orientation: Major axis orientation in clockwise direction as radians. + # expected orientation is calculated frm the minor axis, so convert + orientation = np.abs(np.rad2deg(orientation)) + + if expected_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_max_projection_with_rotation_out_of_plane( + max_projection, + edges, + yc, + xc, + a, + b, + orientation, + f"debug/error_{plane_angle}_{expected_orientation}.png", + ) + assert ( + False + ), f"Major and minor axes should be close to each other, instead got {a} and {b}" + elif expected_orientation is not None: + if not np.allclose(orientation, expected_orientation, atol=atol): + plot_max_projection_with_rotation_out_of_plane( + max_projection, + edges, + yc, + xc, + a, + b, + orientation, + f"debug/error_{plane_angle}_{expected_orientation}.png", + ) + assert ( + False + ), f"Orientation should be close to {expected_orientation}, instead got {orientation}" else: - assert o < 5, "Orientation should be close to 0" + # 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_max_projection_with_rotation_out_of_plane( + max_projection, + edges, + yc, + xc, + a, + b, + orientation, + f"debug/error_{plane_angle}_{expected_orientation}.png", + ) + assert ( + False + ), f"Difference between major and minor axes should be close to {expected_b}, instead got {b}" + + +def plot_max_projection_with_rotation_out_of_plane( + max_projection, edges, yc, xc, a, b, o, title +): + # plot for debugging purposes + + fig, ax = plt.subplots(1, 3, 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[2].imshow(empty_image, cmap="gray") + ax[2].set_title("Ellipse") + ax[2].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 )