Fix vessel deformation bug

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

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)
 

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)
 

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])
+

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])