From 1ce860d1b2d890d8c51128e0e20ea56fef7e34e3 Mon Sep 17 00:00:00 2001 From: faberno Date: Thu, 7 Nov 2024 16:11:42 +0100 Subject: [PATCH] - fixed vessel deformation bug - catched cosine_scaling_factor = 0 case --- simpa/utils/deformation_manager.py | 11 +-- .../CircularTubularStructure.py | 6 +- .../EllipticalTubularStructure.py | 6 +- .../structure_tests/test_elliptical_tubes.py | 34 ++++++++- .../structure_tests/test_tubes.py | 73 +++++++++++++------ 5 files changed, 97 insertions(+), 33 deletions(-) diff --git a/simpa/utils/deformation_manager.py b/simpa/utils/deformation_manager.py index 0c786a0d..4b863026 100644 --- a/simpa/utils/deformation_manager.py +++ b/simpa/utils/deformation_manager.py @@ -25,11 +25,12 @@ def create_deformation_settings(bounds_mm, maximum_z_elevation_mm=1, filter_sigm y_positions_vector = np.linspace(bounds_mm[1][0], bounds_mm[1][1], number_of_boundary_points[1]) # Add random permutations to the y-axis of the division knots - all_scaling_value = np.multiply.outer( - np.cos(x_positions_vector / (bounds_mm[0][1] * (cosine_scaling_factor / - np.pi)) - np.pi / (cosine_scaling_factor * 2)) ** 2, - np.cos(y_positions_vector / (bounds_mm[1][1] * (cosine_scaling_factor / np.pi)) - np.pi / (cosine_scaling_factor * 2)) ** 2) - surface_elevations *= all_scaling_value + if cosine_scaling_factor != 0: + all_scaling_value = np.multiply.outer( + np.cos(x_positions_vector / (bounds_mm[0][1] * (cosine_scaling_factor / + np.pi)) - np.pi / (cosine_scaling_factor * 2)) ** 2, + np.cos(y_positions_vector / (bounds_mm[1][1] * (cosine_scaling_factor / np.pi)) - np.pi / (cosine_scaling_factor * 2)) ** 2) + surface_elevations *= all_scaling_value # This rescales and sets the maximum to 0. surface_elevations = surface_elevations * maximum_z_elevation_mm diff --git a/simpa/utils/libraries/structure_library/CircularTubularStructure.py b/simpa/utils/libraries/structure_library/CircularTubularStructure.py index 3e2fcaa3..ebab63c1 100644 --- a/simpa/utils/libraries/structure_library/CircularTubularStructure.py +++ b/simpa/utils/libraries/structure_library/CircularTubularStructure.py @@ -74,11 +74,11 @@ def get_enclosed_indices(self): torch.arange(self.volume_dimensions_voxels[1]) * self.voxel_spacing, indexing='ij') deformation_values_mm = self.deformation_functional_mm(eval_points) deformation_values_mm = deformation_values_mm.reshape(self.volume_dimensions_voxels[0], - self.volume_dimensions_voxels[1], 1, 1) + self.volume_dimensions_voxels[1], 1) deformation_values_mm = torch.tile(torch.as_tensor( - deformation_values_mm, dtype=torch.float, device=self.torch_device), (1, 1, self.volume_dimensions_voxels[2], 3)) + deformation_values_mm, dtype=torch.float, device=self.torch_device), (1, 1, self.volume_dimensions_voxels[2])) deformation_values_mm /= self.voxel_spacing - target_vector += deformation_values_mm + target_vector[..., -1] += deformation_values_mm del deformation_values_mm cylinder_vector = torch.subtract(end_voxels, start_voxels) diff --git a/simpa/utils/libraries/structure_library/EllipticalTubularStructure.py b/simpa/utils/libraries/structure_library/EllipticalTubularStructure.py index b6b1cf4d..d907dffe 100644 --- a/simpa/utils/libraries/structure_library/EllipticalTubularStructure.py +++ b/simpa/utils/libraries/structure_library/EllipticalTubularStructure.py @@ -81,11 +81,11 @@ def get_enclosed_indices(self): torch.arange(self.volume_dimensions_voxels[1], dtype=torch.float) * self.voxel_spacing, indexing='ij') deformation_values_mm = self.deformation_functional_mm(eval_points) deformation_values_mm = deformation_values_mm.reshape(self.volume_dimensions_voxels[0], - self.volume_dimensions_voxels[1], 1, 1) + self.volume_dimensions_voxels[1], 1) deformation_values_mm = torch.tile(torch.as_tensor( - deformation_values_mm, device=self.torch_device), (1, 1, self.volume_dimensions_voxels[2], 3)) + deformation_values_mm, device=self.torch_device), (1, 1, self.volume_dimensions_voxels[2])) deformation_values_mm /= self.voxel_spacing - target_vector += deformation_values_mm + target_vector[..., -1] += deformation_values_mm del deformation_values_mm cylinder_vector = torch.subtract(end_voxels, start_voxels) diff --git a/simpa_tests/automatic_tests/structure_tests/test_elliptical_tubes.py b/simpa_tests/automatic_tests/structure_tests/test_elliptical_tubes.py index e20681ee..fd77e557 100644 --- a/simpa_tests/automatic_tests/structure_tests/test_elliptical_tubes.py +++ b/simpa_tests/automatic_tests/structure_tests/test_elliptical_tubes.py @@ -5,7 +5,7 @@ import unittest import numpy as np from simpa.utils.libraries.tissue_library import TISSUE_LIBRARY -from simpa.utils import Tags +from simpa.utils import Tags, create_deformation_settings from simpa.utils.settings import Settings from simpa.utils.libraries.structure_library import EllipticalTubularStructure @@ -122,3 +122,35 @@ def test_elliptical_tube_structure_partial_volume_with_eccentricity(self): assert 0 < ets.geometrical_volume[50, 50, 61] < 1 assert 0 < ets.geometrical_volume[30, 50, 50] < 1 assert 0 < ets.geometrical_volume[69, 50, 50] < 1 + + def test_elliptical_tube_structure_deformation(self): + self.global_settings[Tags.DIM_VOLUME_X_MM] = 5 + self.global_settings[Tags.DIM_VOLUME_Y_MM] = 5 + self.global_settings[Tags.DIM_VOLUME_Z_MM] = 5 + self.elliptical_tube_settings[Tags.STRUCTURE_START_MM] = [2.5, 0, 2.5] + self.elliptical_tube_settings[Tags.STRUCTURE_END_MM] = [2.5, 5, 2.5] + self.elliptical_tube_settings[Tags.STRUCTURE_RADIUS_MM] = 0.4 + self.elliptical_tube_settings[Tags.STRUCTURE_ECCENTRICITY] = -0.95 + self.global_settings.set_volume_creation_settings( + { + Tags.STRUCTURES: self.elliptical_tube_settings, + Tags.SIMULATE_DEFORMED_LAYERS: True, + Tags.DEFORMED_LAYERS_SETTINGS: create_deformation_settings( + bounds_mm=[[0, self.global_settings[Tags.DIM_VOLUME_X_MM]], + [0, self.global_settings[Tags.DIM_VOLUME_Y_MM]]], + maximum_z_elevation_mm=2, + filter_sigma=0, + cosine_scaling_factor=0) + } + ) + ets = EllipticalTubularStructure(self.global_settings, self.elliptical_tube_settings) + + # check that vessel was only deformed along z-axis, so there should only be values in the 3 middle slices + assert np.any(ets.geometrical_volume[1:4]) + assert np.all(ets.geometrical_volume[(0, 4), :, :] == 0) + + # check that the vessel was deformed (different from non-deformed case) + center_mask = np.zeros(ets.geometrical_volume.shape, dtype=bool) + center_mask[1:4, :, 2] = True + assert np.any(ets.geometrical_volume[~center_mask]) + diff --git a/simpa_tests/automatic_tests/structure_tests/test_tubes.py b/simpa_tests/automatic_tests/structure_tests/test_tubes.py index 7805dd64..e7c08425 100644 --- a/simpa_tests/automatic_tests/structure_tests/test_tubes.py +++ b/simpa_tests/automatic_tests/structure_tests/test_tubes.py @@ -5,7 +5,7 @@ import unittest import numpy as np from simpa.utils.libraries.tissue_library import TISSUE_LIBRARY -from simpa.utils import Tags +from simpa.utils import Tags, create_deformation_settings from simpa.utils.settings import Settings from simpa.utils.libraries.structure_library import CircularTubularStructure @@ -53,37 +53,40 @@ def test_tube_structures_partial_volume_within_one_voxel(self): self.tube_settings[Tags.STRUCTURE_END_MM] = [0.5, 5, 0.5] self.tube_settings[Tags.STRUCTURE_RADIUS_MM] = 0.4 ts = CircularTubularStructure(self.global_settings, self.tube_settings) - assert 0 < ts.geometrical_volume[0, 0, 0] < 1 - assert 0 < ts.geometrical_volume[0, 4, 0] < 1 + + target_mask = np.zeros(ts.geometrical_volume.shape, dtype=bool) + target_mask[0, :, 0] = True + assert np.all((0 < ts.geometrical_volume[target_mask]) & (ts.geometrical_volume[target_mask] < 1)) + assert np.all(0 == ts.geometrical_volume[~target_mask]) def test_tube_structure_partial_volume_within_two_voxels(self): self.tube_settings[Tags.STRUCTURE_START_MM] = [1, 0, 1] self.tube_settings[Tags.STRUCTURE_END_MM] = [1, 5, 1] self.tube_settings[Tags.STRUCTURE_RADIUS_MM] = 0.4 ts = CircularTubularStructure(self.global_settings, self.tube_settings) - assert 0 < ts.geometrical_volume[0, 0, 0] < 1 - assert 0 < ts.geometrical_volume[1, 0, 0] < 1 - assert 0 < ts.geometrical_volume[0, 1, 0] < 1 - assert 0 < ts.geometrical_volume[0, 0, 1] < 1 - assert 0 < ts.geometrical_volume[1, 1, 0] < 1 - assert 0 < ts.geometrical_volume[1, 0, 1] < 1 - assert 0 < ts.geometrical_volume[0, 1, 1] < 1 - assert 0 < ts.geometrical_volume[1, 1, 1] < 1 + + target_mask = np.zeros(ts.geometrical_volume.shape, dtype=bool) + target_mask[0:2, :, 0:2] = True + assert np.all((0 < ts.geometrical_volume[target_mask]) & (ts.geometrical_volume[target_mask] < 1)) + assert np.all(0 == ts.geometrical_volume[~target_mask]) def test_tube_structure_partial_volume_with_one_full_voxel(self): self.tube_settings[Tags.STRUCTURE_START_MM] = [1.5, 0, 1.5] self.tube_settings[Tags.STRUCTURE_END_MM] = [1.5, 5, 1.5] - self.tube_settings[Tags.STRUCTURE_RADIUS_MM] = np.sqrt(3*0.25) # diagonal of exactly one voxel + self.tube_settings[Tags.STRUCTURE_RADIUS_MM] = np.sqrt(3 * 0.25) # diagonal of exactly one voxel ts = CircularTubularStructure(self.global_settings, self.tube_settings) - assert ts.geometrical_volume[1, 1, 1] == 1 - assert ts.geometrical_volume[0, 1, 0] == 0 - assert ts.geometrical_volume[0, 1, 2] == 0 - assert ts.geometrical_volume[2, 1, 0] == 0 - assert ts.geometrical_volume[2, 2, 2] == 0 - assert 0 < ts.geometrical_volume[0, 1, 1] < 1 - assert 0 < ts.geometrical_volume[1, 1, 0] < 1 - assert 0 < ts.geometrical_volume[1, 1, 2] < 1 - assert 0 < ts.geometrical_volume[2, 1, 1] < 1 + + # check values that should be 1 + assert np.all(ts.geometrical_volume[1, :, 1] == 1) + + # check values that should be fuzzy + target_mask = np.zeros(ts.geometrical_volume.shape, dtype=bool) + target_mask[(0, 1, 1, 2), :, (1, 0, 2, 1)] = True + assert np.all((0 < ts.geometrical_volume[target_mask]) & (ts.geometrical_volume[target_mask] < 1)) + + # check values that should be 0 + target_mask[1, :, 1] = True + assert np.all(0 == ts.geometrical_volume[~target_mask]) def test_tube_structure_partial_volume_across_volume(self): self.tube_settings[Tags.STRUCTURE_START_MM] = [0, 0, 0] @@ -102,3 +105,31 @@ def test_tube_structure_partial_volume_across_volume(self): assert 0 < ts.geometrical_volume[1, 2, 2] < 1 assert 0 < ts.geometrical_volume[1, 2, 3] < 1 assert ts.geometrical_volume[4, 4, 4] == 1 + + def test_tube_structure_deformation(self): + self.tube_settings[Tags.STRUCTURE_START_MM] = [2.5, 0, 2.5] + self.tube_settings[Tags.STRUCTURE_END_MM] = [2.5, 5, 2.5] + self.tube_settings[Tags.STRUCTURE_RADIUS_MM] = 0.4 + self.global_settings.set_volume_creation_settings( + { + Tags.STRUCTURES: self.tube_settings, + Tags.SIMULATE_DEFORMED_LAYERS: True, + Tags.DEFORMED_LAYERS_SETTINGS: create_deformation_settings( + bounds_mm=[[0, self.global_settings[Tags.DIM_VOLUME_X_MM]], + [0, self.global_settings[Tags.DIM_VOLUME_Y_MM]]], + maximum_z_elevation_mm=3, + filter_sigma=0, + cosine_scaling_factor=0) + } + ) + + ts = CircularTubularStructure(self.global_settings, self.tube_settings) + + # check that vessel was only deformed along z-axis, so there should only be values in 1 slice + assert np.any(ts.geometrical_volume[2]) + assert np.all(ts.geometrical_volume[(0, 1, 3, 4), :, :] == 0) + + # check that the vessel was deformed (different from non-deformed case) + center_mask = np.zeros(ts.geometrical_volume.shape, dtype=bool) + center_mask[2, :, 2] = True + assert np.any(ts.geometrical_volume[~center_mask])