Merge remote-tracking branch 'origin/develop'

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "simpa"
-version = "0.8.17"
+version = "0.8.18"
 description = "Simulation and Image Processing for Photonics and Acoustics"
 authors = [
     "Division of Intelligent Medical Systems (IMSY), DKFZ <k.dreher@dkfz-heidelberg.de>",

simpa/core/simulation_modules/volume_creation_module/volume_creation_module_model_based_adapter.py

@@ -7,6 +7,7 @@
 from simpa.utils import Tags
 import numpy as np
 from simpa.utils import create_deformation_settings
+import torch
 
 
 class ModelBasedVolumeCreationAdapter(VolumeCreatorModuleBase):
@@ -99,4 +100,7 @@ def create_simulation_volume(self) -> dict:
 
             global_volume_fractions[mask] += added_volume_fraction[mask]
 
+        # explicitly empty cache to free reserved GPU memory after volume creation
+        torch.cuda.empty_cache()
+
         return volumes

simpa/utils/libraries/structure_library/CircularTubularStructure.py

@@ -1,6 +1,8 @@
 # SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ
 # SPDX-FileCopyrightText: 2021 Janek Groehl
 # SPDX-License-Identifier: MIT
+
+import torch
 import numpy as np
 
 from simpa.utils import Tags
@@ -31,9 +33,9 @@ class CircularTubularStructure(GeometricalStructure):
     """
 
     def get_params_from_settings(self, single_structure_settings):
-        params = (np.asarray(single_structure_settings[Tags.STRUCTURE_START_MM]),
-                  np.asarray(single_structure_settings[Tags.STRUCTURE_END_MM]),
-                  np.asarray(single_structure_settings[Tags.STRUCTURE_RADIUS_MM]),
+        params = (single_structure_settings[Tags.STRUCTURE_START_MM],
+                  single_structure_settings[Tags.STRUCTURE_END_MM],
+                  single_structure_settings[Tags.STRUCTURE_RADIUS_MM],
                   single_structure_settings[Tags.CONSIDER_PARTIAL_VOLUME])
         return params
 
@@ -46,14 +48,17 @@ def to_settings(self):
 
     def get_enclosed_indices(self):
         start_mm, end_mm, radius_mm, partial_volume = self.params
+        start_mm = torch.tensor(start_mm, dtype=torch.float).to(self.torch_device)
+        end_mm = torch.tensor(end_mm, dtype=torch.float).to(self.torch_device)
+        radius_mm = torch.tensor(radius_mm, dtype=torch.float).to(self.torch_device)
         start_voxels = start_mm / self.voxel_spacing
         end_voxels = end_mm / self.voxel_spacing
         radius_voxels = radius_mm / self.voxel_spacing
 
-        x, y, z = np.meshgrid(np.arange(self.volume_dimensions_voxels[0]),
-                              np.arange(self.volume_dimensions_voxels[1]),
-                              np.arange(self.volume_dimensions_voxels[2]),
-                              indexing='ij')
+        x, y, z = torch.meshgrid(torch.arange(self.volume_dimensions_voxels[0]).to(self.torch_device),
+                                 torch.arange(self.volume_dimensions_voxels[1]).to(self.torch_device),
+                                 torch.arange(self.volume_dimensions_voxels[2]).to(self.torch_device),
+                                 indexing='ij')
 
         x = x + 0.5
         y = y + 0.5
@@ -64,30 +69,32 @@ def get_enclosed_indices(self):
         else:
             radius_margin = 0.7071
 
-        target_vector = np.subtract(np.stack([x, y, z], axis=-1), start_voxels)
+        target_vector = torch.subtract(torch.stack([x, y, z], axis=-1), start_voxels)
         if self.do_deformation:
             # the deformation functional needs mm as inputs and returns the result in reverse indexing order...
-            deformation_values_mm = self.deformation_functional_mm(np.arange(self.volume_dimensions_voxels[0], step=1) *
+            deformation_values_mm = self.deformation_functional_mm(torch.arange(self.volume_dimensions_voxels[0]) *
                                                                    self.voxel_spacing,
-                                                                   np.arange(self.volume_dimensions_voxels[1], step=1) *
+                                                                   torch.arange(self.volume_dimensions_voxels[1]) *
                                                                    self.voxel_spacing).T
             deformation_values_mm = deformation_values_mm.reshape(self.volume_dimensions_voxels[0],
                                                                   self.volume_dimensions_voxels[1], 1, 1)
-            deformation_values_mm = np.tile(deformation_values_mm, (1, 1, self.volume_dimensions_voxels[2], 3))
-            target_vector = target_vector + (deformation_values_mm / self.voxel_spacing)
-        cylinder_vector = np.subtract(end_voxels, start_voxels)
+            deformation_values_mm = torch.tile(torch.from_numpy(deformation_values_mm).to(
+                self.torch_device), (1, 1, self.volume_dimensions_voxels[2], 3))
+            target_vector = (target_vector + (deformation_values_mm / self.voxel_spacing)).float()
+        cylinder_vector = torch.subtract(end_voxels, start_voxels)
 
-        target_radius = np.linalg.norm(target_vector, axis=-1) * np.sin(
-            np.arccos((np.dot(target_vector, cylinder_vector)) /
-                      (np.linalg.norm(target_vector, axis=-1) * np.linalg.norm(cylinder_vector))))
+        target_radius = torch.linalg.norm(target_vector, axis=-1) * torch.sin(
+            torch.arccos((torch.matmul(target_vector, cylinder_vector)) /
+                         (torch.linalg.norm(target_vector, axis=-1) * torch.linalg.norm(cylinder_vector))))
 
-        volume_fractions = np.zeros(self.volume_dimensions_voxels)
+        volume_fractions = torch.zeros(tuple(self.volume_dimensions_voxels), dtype=torch.float).to(self.torch_device)
 
         filled_mask = target_radius <= radius_voxels - 1 + radius_margin
         border_mask = (target_radius > radius_voxels - 1 + radius_margin) & \
                       (target_radius < radius_voxels + 2 * radius_margin)
 
         volume_fractions[filled_mask] = 1
+
         volume_fractions[border_mask] = 1 - (target_radius - (radius_voxels - radius_margin))[border_mask]
         volume_fractions[volume_fractions < 0] = 0
 
@@ -96,7 +103,7 @@ def get_enclosed_indices(self):
         else:
             mask = filled_mask
 
-        return mask, volume_fractions[mask]
+        return mask.cpu().numpy(), volume_fractions[mask].cpu().numpy()
 
 
 def define_circular_tubular_structure_settings(tube_start_mm: list,

simpa/utils/libraries/structure_library/EllipticalTubularStructure.py

@@ -2,7 +2,7 @@
 # SPDX-FileCopyrightText: 2021 Janek Groehl
 # SPDX-License-Identifier: MIT
 
-import numpy as np
+import torch
 
 from simpa.utils import Tags
 from simpa.utils.libraries.molecule_library import MolecularComposition
@@ -35,11 +35,11 @@ class EllipticalTubularStructure(GeometricalStructure):
     """
 
     def get_params_from_settings(self, single_structure_settings):
-        params = (np.asarray(single_structure_settings[Tags.STRUCTURE_START_MM]),
-                  np.asarray(single_structure_settings[Tags.STRUCTURE_END_MM]),
-                  np.asarray(single_structure_settings[Tags.STRUCTURE_RADIUS_MM]),
-                  np.asarray(single_structure_settings[Tags.STRUCTURE_ECCENTRICITY]),
-                  np.asarray(single_structure_settings[Tags.CONSIDER_PARTIAL_VOLUME]))
+        params = (single_structure_settings[Tags.STRUCTURE_START_MM],
+                  single_structure_settings[Tags.STRUCTURE_END_MM],
+                  single_structure_settings[Tags.STRUCTURE_RADIUS_MM],
+                  single_structure_settings[Tags.STRUCTURE_ECCENTRICITY],
+                  single_structure_settings[Tags.CONSIDER_PARTIAL_VOLUME])
         return params
 
     def to_settings(self):
@@ -53,15 +53,20 @@ def to_settings(self):
 
     def get_enclosed_indices(self):
         start_mm, end_mm, radius_mm, eccentricity, partial_volume = self.params
+        start_mm = torch.tensor(start_mm, dtype=torch.float).to(self.torch_device)
+        end_mm = torch.tensor(end_mm, dtype=torch.float).to(self.torch_device)
+        radius_mm = torch.tensor(radius_mm, dtype=torch.float).to(self.torch_device)
+        eccentricity = torch.tensor(eccentricity, dtype=torch.float).to(self.torch_device)
+        partial_volume = torch.tensor(partial_volume, dtype=torch.float).to(self.torch_device)
 
         start_voxels = start_mm / self.voxel_spacing
         end_voxels = end_mm / self.voxel_spacing
         radius_voxels = radius_mm / self.voxel_spacing
 
-        x, y, z = np.meshgrid(np.arange(self.volume_dimensions_voxels[0]),
-                              np.arange(self.volume_dimensions_voxels[1]),
-                              np.arange(self.volume_dimensions_voxels[2]),
-                              indexing='ij')
+        x, y, z = torch.meshgrid(torch.arange(self.volume_dimensions_voxels[0]).to(self.torch_device),
+                                 torch.arange(self.volume_dimensions_voxels[1]).to(self.torch_device),
+                                 torch.arange(self.volume_dimensions_voxels[2]).to(self.torch_device),
+                                 indexing='ij')
 
         x = x + 0.5
         y = y + 0.5
@@ -72,40 +77,41 @@ def get_enclosed_indices(self):
         else:
             radius_margin = 0.7071
 
-        target_vector = np.subtract(np.stack([x, y, z], axis=-1), start_voxels)
+        target_vector = torch.subtract(torch.stack([x, y, z], axis=-1), start_voxels)
         if self.do_deformation:
             # the deformation functional needs mm as inputs and returns the result in reverse indexing order...
-            deformation_values_mm = self.deformation_functional_mm(np.arange(self.volume_dimensions_voxels[0], step=1) *
+            deformation_values_mm = self.deformation_functional_mm(torch.arange(self.volume_dimensions_voxels[0]) *
                                                                    self.voxel_spacing,
-                                                                   np.arange(self.volume_dimensions_voxels[1], step=1) *
+                                                                   torch.arange(self.volume_dimensions_voxels[1]) *
                                                                    self.voxel_spacing).T
             deformation_values_mm = deformation_values_mm.reshape(self.volume_dimensions_voxels[0],
                                                                   self.volume_dimensions_voxels[1], 1, 1)
-            deformation_values_mm = np.tile(deformation_values_mm, (1, 1, self.volume_dimensions_voxels[2], 3))
-            target_vector = target_vector + (deformation_values_mm / self.voxel_spacing)
-        cylinder_vector = np.subtract(end_voxels, start_voxels)
+            deformation_values_mm = torch.tile(torch.from_numpy(deformation_values_mm).to(
+                self.torch_device), (1, 1, self.volume_dimensions_voxels[2], 3))
+            target_vector = (target_vector + (deformation_values_mm / self.voxel_spacing)).float()
+        cylinder_vector = torch.subtract(end_voxels, start_voxels)
 
         main_axis_length = radius_voxels/(1-eccentricity**2)**0.25
-        main_axis_vector = np.array([cylinder_vector[1], -cylinder_vector[0], 0])
-        main_axis_vector = main_axis_vector/np.linalg.norm(main_axis_vector) * main_axis_length
+        main_axis_vector = torch.tensor([cylinder_vector[1], -cylinder_vector[0], 0]).to(self.torch_device)
+        main_axis_vector = main_axis_vector/torch.linalg.norm(main_axis_vector) * main_axis_length
 
-        minor_axis_length = main_axis_length*np.sqrt(1-eccentricity**2)
-        minor_axis_vector = np.cross(cylinder_vector, main_axis_vector)
-        minor_axis_vector = minor_axis_vector / np.linalg.norm(minor_axis_vector) * minor_axis_length
+        minor_axis_length = main_axis_length*torch.sqrt(1-eccentricity**2)
+        minor_axis_vector = torch.cross(cylinder_vector, main_axis_vector)
+        minor_axis_vector = minor_axis_vector / torch.linalg.norm(minor_axis_vector) * minor_axis_length
 
-        dot_product = np.dot(target_vector, cylinder_vector)/np.linalg.norm(cylinder_vector)
+        dot_product = torch.matmul(target_vector, cylinder_vector)/torch.linalg.norm(cylinder_vector)
 
-        target_vector_projection = np.multiply(dot_product[:, :, :, np.newaxis], cylinder_vector)
+        target_vector_projection = torch.multiply(dot_product[:, :, :, None], cylinder_vector)
         target_vector_from_projection = target_vector - target_vector_projection
 
-        main_projection = np.dot(target_vector_from_projection, main_axis_vector) / main_axis_length
+        main_projection = torch.matmul(target_vector_from_projection, main_axis_vector) / main_axis_length
 
-        minor_projection = np.dot(target_vector_from_projection, minor_axis_vector) / minor_axis_length
+        minor_projection = torch.matmul(target_vector_from_projection, minor_axis_vector) / minor_axis_length
 
-        radius_crit = np.sqrt(((main_projection/main_axis_length)**2 + (minor_projection/minor_axis_length)**2) *
-                              radius_voxels**2)
+        radius_crit = torch.sqrt(((main_projection/main_axis_length)**2 + (minor_projection/minor_axis_length)**2) *
+                                 radius_voxels**2)
 
-        volume_fractions = np.zeros(self.volume_dimensions_voxels)
+        volume_fractions = torch.zeros(tuple(self.volume_dimensions_voxels), dtype=torch.float).to(self.torch_device)
         filled_mask = radius_crit <= radius_voxels - 1 + radius_margin
         border_mask = (radius_crit > radius_voxels - 1 + radius_margin) & \
                       (radius_crit < radius_voxels + 2 * radius_margin)
@@ -121,7 +127,7 @@ def get_enclosed_indices(self):
         else:
             mask = filled_mask
 
-        return mask, volume_fractions[mask]
+        return mask.cpu().numpy(), volume_fractions[mask].cpu().numpy()
 
 
 def define_elliptical_tubular_structure_settings(tube_start_mm: list,

simpa/utils/libraries/structure_library/HorizontalLayerStructure.py

@@ -2,7 +2,7 @@
 # SPDX-FileCopyrightText: 2021 Janek Groehl
 # SPDX-License-Identifier: MIT
 
-import numpy as np
+import torch
 
 from simpa.utils import Tags
 from simpa.utils.libraries.molecule_library import MolecularComposition
@@ -41,8 +41,8 @@ def to_settings(self):
         return settings
 
     def get_enclosed_indices(self):
-        start_mm = np.asarray(self.params[0])
-        end_mm = np.asarray(self.params[1])
+        start_mm = torch.tensor(self.params[0], dtype=torch.float).to(self.torch_device)
+        end_mm = torch.tensor(self.params[1], dtype=torch.float).to(self.torch_device)
         partial_volume = self.params[2]
         start_voxels = start_mm / self.voxel_spacing
         direction_mm = end_mm - start_mm
@@ -51,32 +51,32 @@ def get_enclosed_indices(self):
         if direction_mm[0] != 0 or direction_mm[1] != 0 or direction_mm[2] == 0:
             raise ValueError("Horizontal Layer structure needs a start and end vector in the form of [0, 0, n].")
 
-        x, y, z = np.meshgrid(np.arange(self.volume_dimensions_voxels[0]),
-                              np.arange(self.volume_dimensions_voxels[1]),
-                              np.arange(self.volume_dimensions_voxels[2]),
-                              indexing='ij')
+        x, y, z = torch.meshgrid(torch.arange(self.volume_dimensions_voxels[0]).to(self.torch_device),
+                                 torch.arange(self.volume_dimensions_voxels[1]).to(self.torch_device),
+                                 torch.arange(self.volume_dimensions_voxels[2]).to(self.torch_device),
+                                 indexing='ij')
 
-        target_vector_voxels = np.subtract(np.stack([x, y, z], axis=-1), start_voxels)
+        target_vector_voxels = torch.subtract(torch.stack([x, y, z], axis=-1), start_voxels)
         target_vector_voxels = target_vector_voxels[:, :, :, 2]
         if self.do_deformation:
             # the deformation functional needs mm as inputs and returns the result in reverse indexing order...
-            deformation_values_mm = self.deformation_functional_mm(np.arange(self.volume_dimensions_voxels[0], step=1) *
+            deformation_values_mm = self.deformation_functional_mm(torch.arange(self.volume_dimensions_voxels[0]) *
                                                                    self.voxel_spacing,
-                                                                   np.arange(self.volume_dimensions_voxels[1], step=1) *
+                                                                   torch.arange(self.volume_dimensions_voxels[1]) *
                                                                    self.voxel_spacing).T
-            target_vector_voxels = target_vector_voxels + (deformation_values_mm.reshape(
+            target_vector_voxels = (target_vector_voxels + torch.from_numpy(deformation_values_mm.reshape(
                 self.volume_dimensions_voxels[0],
-                self.volume_dimensions_voxels[1], 1) / self.voxel_spacing)
+                self.volume_dimensions_voxels[1], 1)).to(self.torch_device) / self.voxel_spacing).float()
 
-        volume_fractions = np.zeros(self.volume_dimensions_voxels)
+        volume_fractions = torch.zeros(tuple(self.volume_dimensions_voxels), dtype=torch.float).to(self.torch_device)
 
         if partial_volume:
             bools_first_layer = ((target_vector_voxels >= -1) & (target_vector_voxels < 0))
 
-            volume_fractions[bools_first_layer] = 1 - np.abs(target_vector_voxels[bools_first_layer])
+            volume_fractions[bools_first_layer] = 1 - torch.abs(target_vector_voxels[bools_first_layer])
 
-            initial_fractions = np.max(volume_fractions, axis=2, keepdims=True)
-            floored_depth_voxels = np.floor(depth_voxels - initial_fractions)
+            initial_fractions = torch.max(volume_fractions, dim=2, keepdims=True)[0]
+            floored_depth_voxels = torch.floor(depth_voxels - initial_fractions)
 
             bools_fully_filled_layers = ((target_vector_voxels >= 0) & (target_vector_voxels < floored_depth_voxels))
 
@@ -96,7 +96,7 @@ def get_enclosed_indices(self):
 
         volume_fractions[bools_fully_filled_layers] = 1
 
-        return bools_all_layers, volume_fractions[bools_all_layers]
+        return bools_all_layers.cpu().numpy(), volume_fractions[bools_all_layers].cpu().numpy()
 
 
 def define_horizontal_layer_structure_settings(molecular_composition: MolecularComposition,